# TODO: Add comment
# 
# Author: rogb
###############################################################################
#Rprof()
rm(list=ls(all=TRUE))
#library(FFT)
library(tseries)
library(RToolbox)



smi <- c("abbn.vx","atln.vx","aden.vx","csgn.vx","holn.vx","baer.vx","lonn.vx","nesn.vx","novn.vx","cfr.vx","rog.vx","sgsn.vx","uhr.vx","slhn.vx","rukn.vx","scmn.vx","synn.vx","syst.vx","ubsn.vx","zurn.vx")
# SMI
# ABB LTD N
symbol <- "abbn.vx"
# ACTELION N
symbol <- "atln.vx"
# ADECCO N (ADEN)
symbol <- "aden.vx"
# CS GROUP N
symbol <- "csgn.vx"
# HOLCIM N
symbol <- "holn.vx"
# JULIUS BAER N
symbol <- "baer.vx"
# LONZA N
symbol <- "lonn.vx"
# NESTLE N
symbol <- "nesn.vx"
# NOVARTIS N
symbol <- "novn.vx"
# RICHEMONT
symbol <- "cfr.vx"
# ROCHE GS
symbol <- "rog.vx"
# SGS N
symbol <- "sgsn.vx"
# SWATCH GROUP I
symbol <- "uhr.vx"
# SWISS LIFE HOLDING N
symbol <- "slhn.vx"
# SWISS RE N
symbol <- "rukn.vx"
# SWISSCOM N
symbol <- "scmn.vx"
# SYNGENTA N
symbol <- "synn.vx"
# SYNTHES N
symbol <- "syst.vx"
# UBS N (UBSN)
symbol <- "ubsn.vx"
# ZURICH FINANCIAL N
symbol <- "zurn.vx"

symbol <- "sbigt1.sw"
symbol <- "sbigt3.sw"
symbol <- "sbigt7.sw"

# Indices
symbol <- "^ssmi" # SMI
symbol <- "^sshi" # SPI
symbol <- "^gspc" # S&P 500
symbol <- "^gdaxi" # DAX
symbol <- "^ftse" # LONDON
symbol <- "^ixic" # NASDAQ Composite
symbol <- "^ndx" # NASDAQ 100
symbol <- "^stoxx50e" # EUROSTOXX 50

symbol <- "^ssmi"

stockCurrent <- read.csv(paste("http://download.finance.yahoo.com/d/quotes.csv?s=",symbol,"&f=sl1d1t1c1ohgv&e=.csv",sep=""),header=F,stringsAsFactors=F)
currentDate <- as.Date(as.character(stockCurrent[3]$V3),format="%m/%d/%Y")
stockData <- get.hist.quote(symbol,quote="AdjClose",compression="m")
stockData[length(stockData)] <- as.numeric(stockCurrent[2])


if(FALSE){
	stockData <- get.hist.quote(symbol,quote="AdjClose",compression="d")
	Y <- year(attributes(stockData)$index)
	M <- month(attributes(stockData)$index)
	futureDates <- getTradingDay(sort(unique(getDateOfNthWeekday(Y,M,5,3))),c("SWISSEXCHANGE","EUREX"),Direction=-1)
	stockData <- stockData[futureDates]
	now <- zoo(as.numeric(stockCurrent[6]),currentDate)
	stockData <- c(stockData,now)
}


# stockData[length(stockData)] <- 105.9

#stockData <- exp(cumsum(c(0,rnorm(131,mean=0-0.08^2/2,sd=0.08))))
# stockData[length(stockData)] <- 268
# stockData <- stockData[-length(stockData)]

# stockData <- exp(cumsum(c(0,arima.sim(list(ma=c(0.5,0.5,0.5,0.5,0.5,0.5)),131)*0.075)))
# ts.plot(stockData)
# sfStop()
now <- proc.time()

stockData.numeric <- as.numeric(stockData)
#stockData.numeric <- c(stockData.numeric,as.numeric(stockCurrent[2]))
#stockData.numeric <- c(stockData.numeric,120)


print(length(stockData.numeric))
print(stockData[length(stockData)])
maMax <- 24
#optimPeriod <- 15
tunePeriod <- 60
nMAused <- 1 # smoothing indicator

optimPeriodMax <- min(c(length(stockData.numeric) - tunePeriod - maMax,round(tunePeriod*1.2)))
#optimPeriodMax <- length(stockData.numeric) - tunePeriod - maMax

optimPeriodMin <- 1
optimPeriodAll <- optimPeriodMin:optimPeriodMax
#optimPeriodAll <- seq(10,100,10)
tuneScore <- c()

for(optimPeriod in optimPeriodAll){
#optimPeriod <- 15
nHist <- optimPeriod + maMax
findMA <- rollApply(stockData.numeric,nHist,function(x,maMax,nMAused){
			library(RToolbox)
			# i <- 0
			# i <- length(stockData.numeric) - nHist
			# x <- stockData.numeric[(1+i):(nHist+i)]
			ma <- 1:maMax
			nMa <- length(ma)-1
			nX <- length(x)
			minLength <- min(nX - ma + 1)
			maData <- sapply(ma,function(maHist){
						tmp <- rollSapply(x,maHist,mean)
						tmp[((nX - maHist + 1)-minLength+1):(nX - maHist + 1)]
					})
			
			lsPos <- t(apply(maData,1,function(x){
								tmp <- rep(1,nMa)
								tmp[x[1]<x[2:(nMa+1)]] <- -1
								tmp[x[1]==x[2:(nMa+1)]] <- 0
								tmp
							}))
			
			stockDiff <- diff(x)[(nX-1-minLength+2):(nX-1)]
			lsPos2 <- lsPos[-minLength,]
			maScore <- as.numeric(t(stockDiff) %*% lsPos2)
			#maScore <- -maScore/apply(stockDiff * lsPos2,2,ES)
			
			names(maScore) <- 2:maMax
				
			maFound <- which(maScore==max(maScore))+1
			
			while(length(maFound)<nMAused){
				maFound <- c(maFound,as.numeric(names(which(maScore[-(maFound-1)]==max(maScore[-(maFound-1)])))))
			}
			
			ls <- lsPos[minLength,maFound-1]
			
			list(ma=maFound,ls=ls,maData=maData[minLength,maFound])
		},maMax,nMAused)

lsPos <- sapply(findMA,function(x)mean(x$ls))
maFound <- sapply(findMA,function(x)mean(x$ma))
maData <- sapply(findMA,function(x)mean(x$maData))

realDiff <- lsPos[-length(findMA)] * diff(stockData.numeric)[-(1:(nHist-1))]
#finalSeries <- (stockData.numeric[nHist] + cumsum(c(0,realDiff))) / stockData.numeric[nHist]

print(paste(optimPeriod,"/",optimPeriodMax,sep=""))
print(sum(realDiff[(length(realDiff)-tunePeriod+1):length(realDiff)]))
#tuneScore <- c(tuneScore,finalSeries[length(finalSeries)]/finalSeries[length(finalSeries)-tunePeriod])
tuneScore <- c(tuneScore,sum(realDiff[(length(realDiff)-tunePeriod+1):length(realDiff)]))
}



optimPeriod1 <- round(mean(optimPeriodAll[which(tuneScore==max(tuneScore))]))
mod <- smooth.spline(optimPeriodAll,tuneScore,df=2*sqrt(length(optimPeriodAll)))
optimPeriod2 <- which(mod$y==max(mod$y))
#optimPeriod2 <- 23
#optimPeriod <- round(mean(c(optimPeriod1,optimPeriod2)))
optimPeriod <- optimPeriod1

nHist <- optimPeriod + maMax
findMA <- rollApply(stockData.numeric,nHist,function(x,maMax,nMAused){
			# i <- 0
			# i <- length(stockData.numeric) - nHist
			# x <- stockData.numeric[(1+i):(nHist+i)]
			ma <- 1:maMax
			nMa <- length(ma)-1
			nX <- length(x)
			minLength <- min(nX - ma + 1)
			maData <- sapply(ma,function(maHist){
						tmp <- rollSapply(x,maHist,mean)
						tmp[((nX - maHist + 1)-minLength+1):(nX - maHist + 1)]
					})
			
			lsPos <- t(apply(maData,1,function(x){
								tmp <- rep(1,nMa)
								tmp[x[1]<x[2:(nMa+1)]] <- -1
								tmp[x[1]==x[2:(nMa+1)]] <- 0
								tmp
							}))
			
			stockDiff <- diff(x)[(nX-1-minLength+2):(nX-1)]
			lsPos2 <- lsPos[-minLength,]
			maScore <- as.numeric(t(stockDiff) %*% lsPos2)
			#maScore <- -maScore/apply(stockDiff * lsPos2,2,ES)
			
			names(maScore) <- 2:maMax
			
			maFound <- which(maScore==max(maScore))+1
			
			while(length(maFound)<nMAused){
				maFound <- c(maFound,as.numeric(names(which(maScore[-(maFound-1)]==max(maScore[-(maFound-1)])))))
			}
			
			ls <- lsPos[minLength,maFound-1]
			
			list(ma=maFound,ls=ls,maData=maData[minLength,maFound])
		},maMax,nMAused)
		
lsPos <- sapply(findMA,function(x)mean(x$ls))
maFound <- sapply(findMA,function(x)mean(x$ma))
maData <- sapply(findMA,function(x)mean(x$maData))

criticalValues <- sapply(findMA[length(findMA)][[1]]$ma,function(ma){
			#ma <- 3
			mean(stockData.numeric[(length(stockData.numeric)-ma+1):(length(stockData.numeric)-1)])
		})

realDiff <- lsPos[-length(findMA)] * diff(stockData.numeric)[-(1:(nHist-1))]
#realDiff <- lsPos * diff(stockData.numeric[-(1:(nHist-2))])

#finalSeries <- (stockData.numeric[nHist] + cumsum(c(0,realDiff))) / stockData.numeric[nHist]
finalSeries <- cumsum(c(0,realDiff))



diffFinal <- diff(finalSeries)
x11()
par(mfrow=c(2,3))

ts.plot(finalSeries)
abline(v=which(diff(lsPos)!=0)+1,lty=2)
abline(v=length(finalSeries)-tunePeriod,col=4,lwd=2)
plot(density(diffFinal))
abline(v=ES(diffFinal),col=2,lwd=3)
abline(v=0)
abline(v=quantile(diffFinal,c(0.25,0.5,0.75)),col=c(3),lty=2)
abline(v=mean(diffFinal),col=2,lwd=3)
text(ES(diffFinal),0,round(ES(diffFinal),3),pos=4)
text(mean(diffFinal),0,round(mean(diffFinal),3),pos=4)

stockDataPlot <- stockData.numeric[(length(stockData.numeric)-length(finalSeries)+1):length(stockData.numeric)]
hist(lsPos,length(unique(lsPos))+1)
#stockDataPlot <- stockDataPlot/stockDataPlot[1]
ts.plot(cbind(stockDataPlot,maData),col=c(1,2),lwd=2,main=symbol)
abline(v=which(diff(lsPos)!=0)+1,lty=2)
abline(v=length(finalSeries)-tunePeriod,col=4,lwd=2)
ts.plot(cbind(maFound,lsPos),col=c(1,2))
abline(v=length(finalSeries)-tunePeriod,col=4,lwd=2)
plot(optimPeriodAll,tuneScore)
lines(mod$x,mod$y)
abline(v=optimPeriod)
text(optimPeriod+1,mean(tuneScore),optimPeriod)

print("Monthly Performance")
print((tuneScore[optimPeriod]/tunePeriod)/stockDataPlot[length(stockDataPlot)]*100)
print(findMA[length(findMA)])
print("currently best ma")
print(maFound[length(maFound)])
print("go long/short")
print(lsPos[length(lsPos)])
print("Optim Period")
print(optimPeriod)
print("critical values")
print(criticalValues)
print("current value")
print(as.numeric(stockCurrent[2]))
print("critical value difference")
print((criticalValues/as.numeric(stockCurrent[2])-1)*100)
print("P/L probability")
print(table(sign(realDiff))/length(realDiff))
print(proc.time()-now)

#summaryRprof()
